Implementacion de modelos de regresion

Supongamos que se quiere predecir para las siguientes 24 horas

import warnings
warnings.filterwarnings("ignore") 

import numpy as np
import pandas as pd
import plotly.express as px
import matplotlib.pyplot as plt
import statsmodels.api as sm
import seaborn as sns

from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.model_selection import TimeSeriesSplit, GridSearchCV
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

from sklearn.neighbors import KNeighborsRegressor
from sklearn.linear_model import Lasso, Ridge
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor
import pickle
from sklearn.svm import SVR


from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy.stats import jarque_bera
from scipy.stats import jarque_bera

data_pm10=pd.read_csv("C:/Users/fonta/OneDrive/Documentos/CDD/6. Sexto semestre/Machine learning/Proyecto final/seriepm10.csv")

data_pm10['pm10_lag1'] = data_pm10['pm10'].shift(1)  
data_pm10['pm10_lag2'] = data_pm10['pm10'].shift(2)  
data_pm10['pm10_lag3'] = data_pm10['pm10'].shift(3)
data_pm10['pm10_lag4'] = data_pm10['pm10'].shift(4)         
data_pm10['pm10_lag5'] = data_pm10['pm10'].shift(5)  
data_pm10 = data_pm10.dropna() 

X = data_pm10[['pm10_lag1', 'pm10_lag2','pm10_lag3','pm10_lag4','pm10_lag5']]
y = data_pm10['pm10']

X.head()

	pm10_lag1	pm10_lag2	pm10_lag3	pm10_lag4	pm10_lag5
5	107.1	56.3	75.6	101.0	61.9
6	92.1	107.1	56.3	75.6	101.0
7	82.5	92.1	107.1	56.3	75.6
8	55.3	82.5	92.1	107.1	56.3
9	72.6	55.3	82.5	92.1	107.1

Metodología: Para cada uno de los modelos de regresión que serán implementados para predecir la particula de PM10 en la estación Ermita se hará una busqueda de los mejores hiperparametros para cada modelo respectivamente. Para probar estos parametros en los conjuntos de test que serán las ultimas obsevaciones de la serie de tiempo, estos conjuntos serán correspondientes de 7, 14, 21 y 28 días con el fin de mostrar como se comporta el modelo en cada una de las ventanas que son de interés.

K-NN

A continuación se realiza un GridSearch con validación cruzada para nuestros datos tomados a traves del tiempo para así hallar el mejor parametro k para nuestro modelo K-NN para predecir la particula de Pm10 en la estación Ermita.

Separando el conjunto de test:

Test={}
ventanas=[7,14,21,28]
for i in ventanas:
    Test[str(i)]= {'trainX': X[:-i*24], 'trainY':y[:-i*24],'xtest':X[-i*24:],'real': y[-i*24:]}

pipeline = Pipeline([
    ('scaler', StandardScaler()),   
    ('knn', KNeighborsRegressor()) 
])

tscv = TimeSeriesSplit(n_splits=5)

param_grid = {
    'knn__n_neighbors': [3,5,7,10,15,20,30]
}

grid_search = GridSearchCV(pipeline, param_grid, cv=tscv, scoring='r2')

grid_search.fit(X, y)

print("Mejores parámetros:", grid_search.best_params_)
print("Score de entrenamineto:", grid_search.best_score_)

#Predicciones en el conjunto de entrenamiento
predicc_knn = grid_search.predict(X)
residuos_knn = y - predicc_knn

Mejores parámetros: {'knn__n_neighbors': 30}
Score de entrenamineto: 0.5648043586722304

Como resultado, obtendremos que el mejor hiperparametro es \(k=30\), el cual es el parámetro con el cual se obtiene el mejor score en el conjunto de entrenamiento. La métrica utilizada fue \(r^2\). Así, para el parámetro \(k=30\) se obtuvo un \(r^2=0.5648%\), sabemos que este score no es el mejor ni el más adecuado, puesto que indica que el modelo k-nn solo alcanza a explicar la variabilidad de los datos en un 56%.

Veamos como se comporta este modelo en los distintos conjuntos de test.

#Conjunto de test
for clave, valor in Test.items():
    xtrain=valor['trainX']
    ytrain=valor['trainY']
    xtest=valor['xtest']
    real=valor['real']
    
    knn=KNeighborsRegressor(30)
    knn.fit(xtrain,ytrain)

    prediccion= knn.predict(xtest)

    Test[clave]['predic']=prediccion
    Test[clave]['residuos']= real-prediccion
    

Errores y métricas

En segundo lugar pasamos a evaluar las métricas, errores y supuestos de nuestro modelo en el conjunto de test, estas nos ayudan a saber un poco mas a detalle como se comporta al momento de predecir los datos, estos errores son medidos por la diferencia entre lo observado y lo predicho, como es un modelo de regresión las métricas evaluadas seran el \(MAE\), \(MSE\), \(RMSE\), \(R2\). Además se evaluaran los supuestos del modelo, siendo estos la independencia y normalidad de los residuos, para esta tarea serán usadas las pruebas analíticas de Ljung-Box y Jarque_bera.

Recordemos las hipotesis para la prueba de independencia Ljung-Box

\[H_0= \textit{Los residuos no estan correlacionados}\]

\[H_1= \textit{Los residuos estan correlacionados}\]

Recordemos las hipotesis para la prueba de independencia Jarque_bera

\[H_0= \textit{Los residuos son normales}\]

\[H_1= \textit{Los residuos no son normales}\]

tabla1 = pd.DataFrame()

for clave, valor in Test.items():
    real = valor['real']
    predic = valor['predic']
    residuos=valor['residuos']
    
    # Calcular las métricas
    mae = mean_absolute_error(real, predic)
    mse = mean_squared_error(real, predic)
    rmse = np.sqrt(mse)
    r2 = r2_score(real, predic)
    Test[clave]['scoretest']=r2

    ljung_box_results = acorr_ljungbox(residuos, lags=[10], return_df=True)
    
    jb_stat, jb_pvalue = jarque_bera(residuos)
    
    tabla_temp = pd.DataFrame({
        'Ventana': [clave],
        'MAE': [mae],
        'MSE': [mse],
        'RMSE': [rmse],
        'R2': [r2],
        'Ljung-Box p-value': [ljung_box_results['lb_pvalue'].values[0]],
        'Jarque-Bera p-value': [jb_pvalue]
    })
    
    tabla1 = pd.concat([tabla1, tabla_temp], ignore_index=True)

tabla1

	Ventana	MAE	MSE	RMSE	R2	Ljung-Box p-value	Jarque-Bera p-value
0	7	6.875582	101.727795	10.086020	0.600924	0.001451	4.873732e-85
1	14	7.349643	111.605169	10.564335	0.621357	0.001527	1.791071e-73
2	21	7.396917	105.644934	10.278372	0.588653	0.000173	1.527648e-65
3	28	7.561096	112.313095	10.597787	0.549661	0.002250	5.740290e-91

En la tabla anterior se muestran los resultados obtenidos en los diferentes grupos de test, si en primer lugar evaluamos \(R^2\) que indica que tanta variabilidad de los datos es captada por el modelo, podemos ver que el mayor \(r^2\) es en el conjunto de 14 días y el menor se encuentra en el de 28 días. Por parte de los errores, podemos apreciar que el MAE mas bajo es tambien en la ventana de 7 días seguido por la de 14 días. Los residuos del modelo en cada una de las ventanas del test no se cumplen, pues bien todos los p-valores son menores que 0.05 lo cual nos lleva a rechazar las \(H_0\).

Veamos si gráficamente se cumplen estas afirmaciones acerca de los residuos.

Gráficos de los supuestos.

sns.set(style="whitegrid")
fig, axes = plt.subplots(4, 2, figsize=(8, 8))  

for i, (k, datos) in enumerate(Test.items()):
    residuo = datos["residuos"]
    ax_hist = axes[i, 0]  
    ax_acf = axes[i, 1]   

    ax_hist.hist(residuo, bins=30, color='#4E148C', edgecolor='black')
    ax_hist.set_title(f"Histograma de los Residuos - ventana {k}")
    ax_hist.set_xlabel("Valor del Residuo")
    ax_hist.set_ylabel("Frecuencia")

    sm.graphics.tsa.plot_acf(residuo, lags=7, ax=ax_acf, color='#4E148C', vlines_kwargs={'color': '#4E148C'}, alpha=0.05)
    ax_acf.set_title(f"Autocorrelación de los Residuos - ventana {k}")

plt.tight_layout()
plt.show()

Analizando la gráfica anterior no se ve muy claramente lo conluído por las pruebas analíticas, pero, si observamos bien en los gráficos de autocorrelación, si hay un pequeño patron entre las barras de la gráfica lo que indica la correlación entre los residuos. Por parte de los histogramas a primera vista se puede apreciar una distribucion parecida a la campana de la normal, pero aun así la prueba analítica afirma que no son normales.

Pasemos a observar gráficamente como son las predicciones realizadas por el modelo en los disitintos conjuntos de test.

Predicciones en el test

sns.set(style="whitegrid")

fig, axes = plt.subplots(4, 1, figsize=(12,10))

for i, (k, datos) in enumerate(Test.items()):
    obs = datos["real"]
    predic = datos["predic"]
    yTrain = datos["trainY"]
    scsoreTest=datos['scoretest']
    ytrain=yTrain[-400:]
    ax = axes[i]
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(obs)), y=obs, label=f"Observado", color='#8e94f2', ax=ax)
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(predic)), y=predic, label=f"Predicho", color="#f77f00", linestyle="--", ax=ax)
    sns.lineplot(x=range(len(ytrain)), y=ytrain, label="Entrenamiento", color="#4E148C", ax=ax)
    
    ax.set_xlabel("")
    ax.set_ylabel("")
    ax.set_title(f"Ventana {k}, score en el test= {round(scsoreTest,2)}")
    ax.legend()

plt.tight_layout()
plt.show()

En las gráficas anteriores pordemos observar como se comportan los distintos conjunto en las prediccioones, sabiendo que el promedio de los score en los conjunto de test esta al rededpr de 60%, podemos ver como estas predicciones siguen la tendencia de los datos reales, no de manera precisa, pero si bastante cerca.

Predicciones en el entrenamiento: Observado vs Predicho

# Crear la figura
fig = px.line(data_pm10, x='fecha', y=[y, predicc_knn], 
              labels={'value': 'PM10', 'variable': 'Leyenda'})

# Asignar nombres a las trazas y actualizar el color y el estilo
fig.data[0].name = 'pm10_observado'  # Serie observada
fig.data[1].name = 'pm10_predicho'  # Serie predicha

# Cambiar el estilo de las líneas
fig.update_traces(line=dict(width=2))  
fig.update_traces(line_color='#9d4edd', selector=dict(name='pm10_observado'))  
fig.update_traces(line_color='#3c096c', selector=dict(name='pm10_predicho')) 

# Configurar el diseño
fig.update_layout(
    title='Predicción vs Observación de PM10: K-NN',
    xaxis_title='Fecha',
    yaxis_title='Concentración de PM10',
    legend_title='Leyenda',
    template='plotly_white',
    hovermode='x unified',
    width=1300,
    height=500,
    legend=dict(x=1, y=1, xanchor='center', yanchor='bottom')
)

# Mostrar la gráfica
fig.show()

De manera mas general, tenemos las predicciones en el conjunto de entrenamiento, ya una vez entrenado el modelo podemos ver como el modelo lograría predecir esos datos de entrenamiento generando como resultado la gráfica anterios, recordando que el score obtenido en el conjunto de entrenamiento fue de \(R^2= 56%\) lo cual indica que el modelo capta de manera moderada la variabilidad de los datos, así como se ve en la gráfica, vemos que al modelo le cuesta predecir los picos del pm10, por ejemplo al rededor de finales de 2019.

Lasso

La misma lógica que se aplico en el modelo anterior será usada para la regresión Lasso, para así hallar el mejor parametro de reguralización \(\alpha\), ademas veremos como se comporta el modelo en los diferentes conjuntos de test.

Separamos nuevamente los diferentes conjuntos de test, recordamos que son los ultimos días, sean 7,14,21,28.

Test_lasso={}
ventanas=[7,14,21,28]
for i in ventanas:
    Test_lasso[str(i)]= {'trainX': X[:-i*24], 'trainY':y[:-i*24],'xtest':X[-i*24:],'real': y[-i*24:]}

A continuación se realiza la hiperparametrización para le modelo Lasso.

pipeline = Pipeline([
    ('scaler', StandardScaler()),   
    ('lasso', Lasso()) 
])
tscv = TimeSeriesSplit(n_splits=5)

param_grid = {
    'lasso__alpha': [0.001, 0.01, 0.1, 1.0, 10.0] 
}

grid_search = GridSearchCV(pipeline, param_grid, cv=tscv, scoring='r2')

grid_search.fit(X, y)

print("Mejores parámetros:", grid_search.best_params_)
print("Score de entrenamiento:", grid_search.best_score_)

#Predicciones en el conjunto de entrenamiento
predicc_lasso = grid_search.predict(X)
residuos_lasso = y - predicc_lasso

Mejores parámetros: {'lasso__alpha': 0.01}
Score de entrenamiento: 0.5649782914088946

Obtuvimos como resultado que de todos los posibles parametros de \(\alpha\) el que nos entrega el mejor score siendo este de 56.49% es el \(\alpha = 0.01\). Podemos ver que este score en el conjunto de entrenamiento no es muy distante del que se obtuvo con el modelo anterior.

Pasemos a evaluar este modelo en los conjuntos de test.

for clave, valor in Test_lasso.items():
    xtrain=valor['trainX']
    ytrain=valor['trainY']
    xtest=valor['xtest']
    real=valor['real']
    
    lasso=Lasso(alpha=0.01)
    lasso.fit(xtrain,ytrain)

    prediccion= lasso.predict(xtest)

    Test_lasso[clave]['predic']=prediccion
    Test_lasso[clave]['residuos']= real-prediccion

Evaluemos los rendimiento en cada conjunto de test mediante las diferentes métricas y supuesto para los residuos al igual que lo hicimos con el modelo anterior.

Errores y métricas

tabla2 = pd.DataFrame()

for clave, valor in Test_lasso.items():
    real = valor['real']
    predic = valor['predic']
    residuos=valor['residuos']
    
    # Calcular las métricas
    mae = mean_absolute_error(real, predic)
    mse = mean_squared_error(real, predic)
    rmse = np.sqrt(mse)
    r2 = r2_score(real, predic)
    Test_lasso[clave]['scoretest']=r2

    ljung_box_results = acorr_ljungbox(residuos, lags=[10], return_df=True)
    
    jb_stat, jb_pvalue = jarque_bera(residuos)
    
    tabla_temp = pd.DataFrame({
        'Ventana': [clave],
        'MAE': [mae],
        'MSE': [mse],
        'RMSE': [rmse],
        'R2': [r2],
        'Ljung-Box p-value': [ljung_box_results['lb_pvalue'].values[0]],
        'Jarque-Bera p-value': [jb_pvalue]
    })
    
    tabla2 = pd.concat([tabla2, tabla_temp], ignore_index=True)

# Mostrar el DataFrame final
tabla2

	Ventana	MAE	MSE	RMSE	R2	Ljung-Box p-value	Jarque-Bera p-value
0	7	7.073974	104.603971	10.227608	0.589641	0.007450	1.859167e-96
1	14	7.448370	113.226365	10.640788	0.615856	0.001458	2.486375e-86
2	21	7.411331	107.403201	10.363552	0.581807	0.000042	4.706713e-81
3	28	7.623725	115.302152	10.737884	0.537676	0.005496	6.055035e-114

Si observamos la tabla con los resultados para cada uno de las ventanas en el conjunto de test, en primer lugar tenemos que el mejor \(R^2\) es de 0.61, este fue obtenido en el conjunto de test de 14 días, en este mismo conjunto se obtuvo un MAE de 7.4 y por parte del RMSE que es la métrica encargada de medir el error en las unidades reales de la variable, en este conjunto se obtuvo un RMSE= 10.64 en las demas ventanas el score estan al rededor de 58%. Ademas por parte de los residuos vemos que en nunguno de los conjutos de test se cumplen los supuestos de normlidad ni de independencia.

A continuación veamos las gráficas de los residuos del modelo en cada uuno de los conjuntos de test.

Gráficos de los supuestos

sns.set(style="whitegrid")
fig, axes = plt.subplots(4, 2, figsize=(8, 8))  

for i, (k, datos) in enumerate(Test_lasso.items()):
    residuo = datos["residuos"]
    ax_hist = axes[i, 0]  
    ax_acf = axes[i, 1]   

    ax_hist.hist(residuo, bins=30, color='#4E148C', edgecolor='black')
    ax_hist.set_title(f"Histograma de los Residuos - ventana {k}")
    ax_hist.set_xlabel("Valor del Residuo")
    ax_hist.set_ylabel("Frecuencia")

    sm.graphics.tsa.plot_acf(residuo, lags=7, ax=ax_acf, color='#4E148C', vlines_kwargs={'color': '#4E148C'}, alpha=0.05)
    ax_acf.set_title(f"Autocorrelación de los Residuos - ventana {k}")

plt.tight_layout()
plt.show()

Las conclusiones que podemos encontrar en las gráficas anteriores es que parece ser que los supuestos de los modelos si se cumplen, pero sabemos que en estadística las pruebas analíticas son las que tienen la ultima palabra, asi que ainque parezca que los residuos tienen una distribución normal y que las correlaciones no son tan significativas, pero ya vimos que las pruebas analíticas confirman que no se cumplen los residuos.

Predicciones en el test

sns.set(style="whitegrid")

fig, axes = plt.subplots(4, 1, figsize=(12,10))

for i, (k, datos) in enumerate(Test_lasso.items()):
    obs = datos["real"]
    predic = datos["predic"]
    yTrain = datos["trainY"]
    scsoreTest=datos['scoretest']
    ytrain=yTrain[-400:]
    ax = axes[i]
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(obs)), y=obs, label=f"Observado", color='#8e94f2', ax=ax)
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(predic)), y=predic, label=f"Predicho", color="#f77f00", linestyle="--", ax=ax)
    sns.lineplot(x=range(len(ytrain)), y=ytrain, label="Entrenamiento", color="#4E148C", ax=ax)
    
    ax.set_xlabel("")
    ax.set_ylabel("")
    ax.set_title(f"Ventana {k}, score en el test= {round(scsoreTest,2)}")
    ax.legend()

plt.tight_layout()
plt.show()

Las gráficas anteriores son utiles para mostrar como se traducen esos scores de los que hemos estado hablando a una forma donde se puede visualizar, en cada conjunto de test el modelo busca seguir de cerca los valores reales, pero vemos que las predicciones fallan un poco en los puntos extremos de la serie, por eso es que vemos que la variabilidad que es captada por el modelo en cada conjunto es moderada.

Predicciones en el entrenamiento: Observado VS Predicho

# Crear la figura
fig = px.line(data_pm10, x='fecha', y=[y, predicc_lasso], 
              labels={'value': 'PM10', 'variable': 'Leyenda'})

# Asignar nombres a las trazas y actualizar el color y el estilo
fig.data[0].name = 'pm10_observado'  # Serie observada
fig.data[1].name = 'pm10_predicho'  # Serie predicha

# Cambiar el estilo de las líneas
fig.update_traces(line=dict(width=2))  
fig.update_traces(line_color='#9d4edd', selector=dict(name='pm10_observado'))  
fig.update_traces(line_color='#3c096c', selector=dict(name='pm10_predicho')) 

# Configurar el diseño
fig.update_layout(
    title='Predicción vs Observación de PM10: Lasso',
    xaxis_title='Fecha',
    yaxis_title='Concentración de PM10',
    legend_title='Leyenda',
    template='plotly_white',
    hovermode='x unified',
    width=1300,
    height=500,
    legend=dict(x=1, y=1, xanchor='center', yanchor='bottom')
)

# Mostrar la gráfica
fig.show()

Si estudiamos como se ven las predicciones de todo el conjunto de entrenamiento, vemos que se cumple lo que mencionabamos anteriormente, el modelo tiene problemas al momento de de predecir los valores extremos, aunque si compaaramos.

Ridge

Siguiendo la metodología, pasamos a implementar el modelo de regresión Ridge y a analizar sus resultados, veamos.

Separamos cada conjunto de test.

Test_ridge={}
ventanas=[7,14,21,28]
for i in ventanas:
    Test_ridge[str(i)]= {'trainX': X[:-i*24], 'trainY':y[:-i*24],'xtest':X[-i*24:],'real': y[-i*24:]}

Hacemos la hiperparametrización respectiva para el modelo Ridge.

pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('ridge', Ridge()) 
])
tscv = TimeSeriesSplit(n_splits=5)

param_grid = {
    'ridge__alpha': [0.001, 0.01, 0.1, 1.0, 10.0, 100.0]  
}

grid_search = GridSearchCV(pipeline, param_grid, cv=tscv, scoring='r2')

grid_search.fit(X, y)

print("Mejores parámetros:", grid_search.best_params_)
print("Score de entrenamiento:", grid_search.best_score_)

predicc_ridge = grid_search.predict(X)

# Calcular los residuos
residuos_ridge = y - predicc_ridge

Mejores parámetros: {'ridge__alpha': 0.001}
Score de entrenamiento: 0.5649771137261234

En este caso para la regresion Ridge encontramos que el parametro de regularización que entrega le mejor score en el conjunto de entrenamiento es \(\alpha = 0.001\), el score en el entrenamiento es de 56.49% que si recordamos es igual al de la regresion Lasso.

Veamos el rendimiento en los conjuntos de test mediante los errores y las métricas.

for clave, valor in Test_ridge.items():
    xtrain=valor['trainX']
    ytrain=valor['trainY']
    xtest=valor['xtest']
    real=valor['real']
    
    ridge=Ridge(alpha=0.001)
    ridge.fit(xtrain,ytrain)

    prediccion= ridge.predict(xtest)

    Test_ridge[clave]['predic']=prediccion
    Test_ridge[clave]['residuos']= real-prediccion

Errores y métricas

tabla3 = pd.DataFrame()

for clave, valor in Test_ridge.items():
    real = valor['real']
    predic = valor['predic']
    residuos=valor['residuos']
    
    # Calcular las métricas
    mae = mean_absolute_error(real, predic)
    mse = mean_squared_error(real, predic)
    rmse = np.sqrt(mse)
    r2 = r2_score(real, predic)
    Test_ridge[clave]['scoretest']=r2

    ljung_box_results = acorr_ljungbox(residuos, lags=[10], return_df=True)
    
    jb_stat, jb_pvalue = jarque_bera(residuos)
    
    tabla_temp = pd.DataFrame({
        'Ventana': [clave],
        'MAE': [mae],
        'MSE': [mse],
        'RMSE': [rmse],
        'R2': [r2],
        'Ljung-Box p-value': [ljung_box_results['lb_pvalue'].values[0]],
        'Jarque-Bera p-value': [jb_pvalue]
    })
    
    tabla3 = pd.concat([tabla3, tabla_temp], ignore_index=True)

# Mostrar el DataFrame final
tabla3

	Ventana	MAE	MSE	RMSE	R2	Ljung-Box p-value	Jarque-Bera p-value
0	7	7.073715	104.599215	10.227376	0.589659	0.007471	1.731642e-96
1	14	7.448146	113.221035	10.640537	0.615874	0.001467	2.411876e-86
2	21	7.411178	107.399532	10.363375	0.581821	0.000043	4.643038e-81
3	28	7.623578	115.299627	10.737766	0.537686	0.005525	5.813875e-114

Estudiando los resultados obtenidos para el modelo Ridge vemos que son bastante similiares a los resultados del modelo Lasso, el mejor score dentro de todos los conjuntos de test es nuevamente es el de los ultimos 14 días con un score de 61% y los errores relativamente bajos y nuevamente los supuestos no se cumplen en este modelo, los residuos en cada conjunto no cumplen con la normalidad y la independencia.

Gráficos de los supuestos

sns.set(style="whitegrid")
fig, axes = plt.subplots(4, 2, figsize=(8, 8))  

for i, (k, datos) in enumerate(Test_ridge.items()):
    residuo = datos["residuos"]
    ax_hist = axes[i, 0]  
    ax_acf = axes[i, 1]   

    ax_hist.hist(residuo, bins=30, color='#4E148C', edgecolor='black')
    ax_hist.set_title(f"Histograma de los Residuos - ventana {k}")
    ax_hist.set_xlabel("Valor del Residuo")
    ax_hist.set_ylabel("Frecuencia")

    sm.graphics.tsa.plot_acf(residuo, lags=7, ax=ax_acf, color='#4E148C', vlines_kwargs={'color': '#4E148C'}, alpha=0.05)
    ax_acf.set_title(f"Autocorrelación de los Residuos - ventana {k}")

plt.tight_layout()
plt.show()

Analizando las gráficas de autocorrelación se nota un ligero patron entre cada uno de los lags de la figura, lo que evidencia la correlación que esta presente en los residuos, además, por parte de los histogramas que parecen mostrar una distribución cercana a la normal pero tenemos de respaldo la prueba analítica de Jarque-Bera que indica que los residuos no son normales.

Veamos las predicciones en cada conjunto de test.

Predicciones en el test

sns.set(style="whitegrid")

fig, axes = plt.subplots(4, 1, figsize=(10,8))

for i, (k, datos) in enumerate(Test_ridge.items()):
    obs = datos["real"]
    predic = datos["predic"]
    yTrain = datos["trainY"]
    scsoreTest=datos['scoretest']
    ytrain=yTrain[-400:]
    ax = axes[i]

    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(obs)), y=obs, label=f"Observado", color='#8e94f2', ax=ax)
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(predic)), y=predic, label=f"Predicho", color="#f77f00", linestyle="--", ax=ax)
    sns.lineplot(x=range(len(ytrain)), y=ytrain, label="Entrenamiento", color="#4E148C", ax=ax)
    
    ax.set_xlabel("")
    ax.set_ylabel("")
    ax.set_title(f"Ventana {k}, score en el test= {round(scsoreTest,2)}")
    ax.legend()

plt.tight_layout()
plt.show()

En la gráfica anterior se tienen las predicciones hechas por el modelo Ridge en cada uno de los conjuntos de test. Vemos que hay entre el score obtenido y que tan bien el modelo logra predecir los extremos de la serie de tiempo, si vemos el score mas bajo se dió en el conjunto de 28 días y la prediccion en los extremos no es tan precisa a comparacion de la ventana de 14 que sabemos es la que cuenta con el mejor score en el test y cuenta con un poco mas de precision en los extremos.

Gráfica predicciones en el entrenamiento: Observado vs Predicho

# Crear la figura
fig = px.line(data_pm10, x='fecha', y=[y, predicc_ridge], 
              labels={'value': 'PM10', 'variable': 'Leyenda'})

# Asignar nombres a las trazas y actualizar el color y el estilo
fig.data[0].name = 'pm10_observado'  # Serie observada
fig.data[1].name = 'pm10_predicho'  # Serie predicha

# Cambiar el estilo de las líneas
fig.update_traces(line=dict(width=2))  
fig.update_traces(line_color='#9d4edd', selector=dict(name='pm10_observado'))  
fig.update_traces(line_color='#3c096c', selector=dict(name='pm10_predicho')) 

# Configurar el diseño
fig.update_layout(
    title='Predicción vs Observación de PM10: Ridge',
    xaxis_title='Fecha',
    yaxis_title='Concentración de PM10',
    legend_title='Leyenda',
    template='plotly_white',
    hovermode='x unified',
    width=1300,
    height=500,
    legend=dict(x=1, y=1, xanchor='center', yanchor='bottom')
)

# Mostrar la gráfica
fig.show()

Como lo mencionamos y era de esperarse al predecir los valores de entrenmiento se presenta el mismo comportamiento, este modelo Ridge tambien esta teniendo problemas al predecir los valores extremos en la serie de tiempo de PM10 en la estacion Ermita.

Random Forest

Veamos el rendimiento del modelo de random forest al momento de predecir la particula del PM10.

Separemos los distintos conjuntos de test.

Test_rf={}
ventanas=[7,14,21,28]
for i in ventanas:
    Test_rf[str(i)]= {'trainX': X[:-i*24], 'trainY':y[:-i*24],'xtest':X[-i*24:],'real': y[-i*24:]}

Encontremos los mejores parametros para este modelo de Random Forest.

pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('random_forest', RandomForestRegressor())
])

tscv = TimeSeriesSplit(n_splits=5)

param_grid = {
    'random_forest__n_estimators': [50, 100],  
    'random_forest__max_depth':  [10, 20],     
    'random_forest__min_samples_split': [2, 5, 10]  
}

# Configurar el GridSearchCV
grid_search = GridSearchCV(pipeline, param_grid, cv=tscv, scoring='r2')

# Ajustar el modelo
grid_search.fit(X, y)

# Imprimir los mejores parámetros y el score de entrenamiento
print("Mejores parámetros:", grid_search.best_params_)
print("Score de entrenamiento:", grid_search.best_score_)

# Realizar predicciones con el mejor modelo
predicc_rf = grid_search.predict(X)

# Calcular los residuos
residuos_rf = y - predicc_rf

with open('modelo_rf.pkl', 'wb') as f:
    pickle.dump({
        'modelo': grid_search.best_estimator_,
        'predicciones': predicc_rf,
        'residuos': residuos_rf
    }, f)

print("Modelo, predicciones y residuos guardados en 'modelo_rf.pkl'")

Mejores parámetros: {'random_forest__max_depth': 10, 'random_forest__min_samples_split': 10, 'random_forest__n_estimators': 50}
Score de entrenamiento: 0.5675789447397046
Modelo, predicciones y residuos guardados en 'modelo_rf.pkl'

Vemos que los mejores parametros y los cuales entregan el mejor score \(R^2= 56%\) son una profundidad maxima para cada arbol de 10, ademas de un minimo en muestra de 10 arboles y un total de 50 estimadores, la combinación de estos parametros entregan el mejor score en el conjunto de train.

with open('modelo_rf.pkl', 'rb') as f:
    model_rf = pickle.load(f)

Veamos como se comporta el modelo en los conjuntos de test estudiando sus métricas y residuos.

for clave, valor in Test_rf.items():
    xtrain = valor['trainX']
    ytrain = valor['trainY']
    xtest = valor['xtest']
    real = valor['real']
    
    rf_model = RandomForestRegressor(
        n_estimators=100,     
        max_depth=10,         
        min_samples_split=10,  
    )
    
    rf_model.fit(xtrain, ytrain)
    
    prediccion = rf_model.predict(xtest)

    Test_rf[clave]['predic'] = prediccion
    Test_rf[clave]['residuos'] = real - prediccion

Métricas y errores

tabla3 = pd.DataFrame()

for clave, valor in Test_rf.items():
    real = valor['real']
    predic = valor['predic']
    residuos=valor['residuos']
    
    # Calcular las métricas
    mae = mean_absolute_error(real, predic)
    mse = mean_squared_error(real, predic)
    rmse = np.sqrt(mse)
    r2 = r2_score(real, predic)
    Test_rf[clave]['scoretest']=r2

    ljung_box_results = acorr_ljungbox(residuos, lags=[10], return_df=True)
    
    jb_stat, jb_pvalue = jarque_bera(residuos)
    
    tabla_temp = pd.DataFrame({
        'Ventana': [clave],
        'MAE': [mae],
        'MSE': [mse],
        'RMSE': [rmse],
        'R2': [r2],
        'Ljung-Box p-value': [ljung_box_results['lb_pvalue'].values[0]],
        'Jarque-Bera p-value': [jb_pvalue]
    })
    
    tabla3 = pd.concat([tabla3, tabla_temp], ignore_index=True)

tabla3

	Ventana	MAE	MSE	RMSE	R2	Ljung-Box p-value	Jarque-Bera p-value
0	7	6.654849	95.545911	9.774759	0.625175	0.013958	1.923931e-153
1	14	7.077893	103.854556	10.190906	0.647652	0.022222	1.162537e-106
2	21	7.158516	100.713200	10.035597	0.607855	0.002082	7.042026e-79
3	28	7.345389	108.613025	10.421757	0.564497	0.009254	2.431053e-118

Estudiando la tabla anterior en la que se muestran los resultados en los conjuntos de test vemos que el mejor score es de 64% y fue obtenida en el conjunto de los ultimos 14 días, donde ademas los errores no son tan grandes al igual que en los demas conjutos. De la misma manera este modelo no cumple con los supuestos de independencia y normalidad en los residuos.

Gráfica para los supuestos

sns.set(style="whitegrid")
fig, axes = plt.subplots(4, 2, figsize=(8, 8))  

for i, (k, datos) in enumerate(Test_rf.items()):
    residuo = datos["residuos"]
    ax_hist = axes[i, 0]  
    ax_acf = axes[i, 1]   

    ax_hist.hist(residuo, bins=30, color='#4E148C', edgecolor='black')
    ax_hist.set_title(f"Histograma de los Residuos - ventana {k}")
    ax_hist.set_xlabel("Valor del Residuo")
    ax_hist.set_ylabel("Frecuencia")

    sm.graphics.tsa.plot_acf(residuo, lags=7, ax=ax_acf, color='#4E148C', vlines_kwargs={'color': '#4E148C'}, alpha=0.05)
    ax_acf.set_title(f"Autocorrelación de los Residuos - ventana {k}")

plt.tight_layout()
plt.show()

En las gráficas anteriores podemos ver los histogramas y autocorrelaciones para los residuos del modelo, aunque parecen seguir una distribución normal, sabemos que la prueba analítica arrojó que no se cumple el supuesto de normalidad y en las autocorrelaciones se puede ver un ligero patrón lo que evidencia la correlacion que existe entre ellos.

Predicciones en el test

sns.set(style="whitegrid")

fig, axes = plt.subplots(4, 1, figsize=(12,8))

for i, (k, datos) in enumerate(Test_rf.items()):
    obs = datos["real"]
    predic = datos["predic"]
    yTrain = datos["trainY"]
    scsoreTest=datos['scoretest']
    ytrain=yTrain[-400:]
    ax = axes[i]
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(obs)), y=obs, label=f"Observado", color='#8e94f2', ax=ax)
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(predic)), y=predic, label=f"Predicho", color="#f77f00", linestyle="--", ax=ax)
    sns.lineplot(x=range(len(ytrain)), y=ytrain, label="Entrenamiento", color="#4E148C", ax=ax)
    
    ax.set_xlabel("")
    ax.set_ylabel("")
    ax.set_title(f"Ventana {k}, score en el test= {round(scsoreTest,2)}")
    ax.legend()

plt.tight_layout()
plt.show()

En la figura anterior podemos ver cada uno de las predicciones hechas en cada uno de los conjuntos de test, vemos que estas predicciones no son muy diferentes a las demas que hemos estudiado en los demas modelos, a pesar que los scores no son tan altos las predicciones siguen de cerca los datos reales.

Predicciones en el conjunto de entrenamiento

# Crear la figura
fig = px.line(data_pm10, x='fecha', y=[y, model_rf['predicciones']], 
              labels={'value': 'PM10', 'variable': 'Leyenda'})

# Asignar nombres a las trazas y actualizar el color y el estilo
fig.data[0].name = 'pm10_observado'  # Serie observada
fig.data[1].name = 'pm10_predicho'  # Serie predicha

# Cambiar el estilo de las líneas
fig.update_traces(line=dict(width=2))  
fig.update_traces(line_color='#9d4edd', selector=dict(name='pm10_observado'))  
fig.update_traces(line_color='#3c096c', selector=dict(name='pm10_predicho')) 


# Configurar el diseño
fig.update_layout(
    title='Predicción vs Observación de PM10: Random Forest',
    xaxis_title='Fecha',
    yaxis_title='Concentración de PM10',
    legend_title='Leyenda',
    template='plotly_white',
    hovermode='x unified',
    width=1000,
    height=500,
    legend=dict(x=1, y=1, xanchor='center', yanchor='bottom')
)

fig.show()

XGBoost

Veamos como es el comportamiento del modelo XGboost para nuestros datos.

Test_xgboost={}
ventanas=[7,14,21,28]
for i in ventanas:
    Test_xgboost[str(i)]= {'trainX': X[:-i*24], 'trainY':y[:-i*24],'xtest':X[-i*24:],'real': y[-i*24:]}

Realizamos la respectiva hiperparametrización para el modelo.

pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('xgboost', XGBRegressor(tree_method='gpu_hist', gpu_id=0))  
])

tscv = TimeSeriesSplit(n_splits=5)

param_grid = {
    'xgboost__n_estimators': [50, 100],      
    'xgboost__max_depth': [3, 6, 10],             
    'xgboost__learning_rate': [0.01, 0.1, 0.2],   
    'xgboost__subsample': [0.8, 1.0]              
}

grid_search = GridSearchCV(pipeline, param_grid, cv=tscv, scoring='r2')


grid_search.fit(X, y)

# Imprimir los mejores parámetros y el score de entrenamiento
print("Mejores parámetros:", grid_search.best_params_)
print("Score de entrenamiento:", grid_search.best_score_)

# Realizar predicciones con el mejor modelo
predicc_xgb = grid_search.predict(X)

# Calcular los residuos
residuos_xgb = y - predicc_xgb

Mejores parámetros: {'xgboost__learning_rate': 0.1, 'xgboost__max_depth': 3, 'xgboost__n_estimators': 50, 'xgboost__subsample': 0.8}
Score de entrenamiento: 0.5731079997865662

El mejor resultado en el conjunto entrenamiento fue obtenido por una tasa de aprendizaje del 0.1, una profundidad maxima de 3, 50 estimadores, y un score en el conjunto de entrenamiento de 57,3%. Veamos el rendimiento en los conjuntos de test.

#conjunto de test
for clave, valor in Test_xgboost.items():
    xtrain = valor['trainX']
    ytrain = valor['trainY']
    xtest = valor['xtest']
    real = valor['real']
    
    xgb_model = XGBRegressor(
        tree_method='gpu_hist',  
        gpu_id=0,
        learning_rate=0.1,
        max_depth=3,
        n_estimators=50,
        subsample=0.8
    )
    
    # Ajustar el modelo
    xgb_model.fit(xtrain, ytrain)
    
    # Realizar la predicción
    prediccion = xgb_model.predict(xtest)
    
    Test_xgboost[clave]['predic'] = prediccion
    Test_xgboost[clave]['residuos'] = real - prediccion

Observemos los errores y métricas para los residuos de este modelo.

Errores y métricas

tabla4 = pd.DataFrame()

for clave, valor in Test_xgboost.items():
    real = valor['real']
    predic = valor['predic']
    residuos=valor['residuos']
    
    # Calcular las métricas
    mae = mean_absolute_error(real, predic)
    mse = mean_squared_error(real, predic)
    rmse = np.sqrt(mse)
    r2 = r2_score(real, predic)
    Test_xgboost[clave]['scoretest']=r2

    ljung_box_results = acorr_ljungbox(residuos, lags=[10], return_df=True)
    
    jb_stat, jb_pvalue = jarque_bera(residuos)
    
    tabla_temp = pd.DataFrame({
        'Ventana': [clave],
        'MAE': [mae],
        'MSE': [mse],
        'RMSE': [rmse],
        'R2': [r2],
        'Ljung-Box p-value': [ljung_box_results['lb_pvalue'].values[0]],
        'Jarque-Bera p-value': [jb_pvalue]
    })
    
    tabla4 = pd.concat([tabla4, tabla_temp], ignore_index=True)

tabla4

	Ventana	MAE	MSE	RMSE	R2	Ljung-Box p-value	Jarque-Bera p-value
0	7	6.862827	96.207918	9.808564	0.622578	0.005394	5.482729e-82
1	14	7.225789	104.512607	10.223141	0.645420	0.002161	1.445462e-68
2	21	7.288449	102.131383	10.106007	0.602333	0.000039	1.306398e-59
3	28	7.491330	110.366223	10.505533	0.557467	0.001509	1.357846e-94

Analizando la tabla obtenida podemos concluir que el mejor obtenido fue para los 14 días con un score de 64.5% y en este caso el menor score fue para los 28 días. Podemos observar que los errores son bastante pequeños aunque el ajuste del modelo no sea el mejor. De nuevo, no se cumplen los supuestos de normalidad e independencia.

Gráfica para los supuestos

sns.set(style="whitegrid")
fig, axes = plt.subplots(4, 2, figsize=(8, 8))  

for i, (k, datos) in enumerate(Test_xgboost.items()):
    residuo = datos["residuos"]
    ax_hist = axes[i, 0]  
    ax_acf = axes[i, 1]   

    ax_hist.hist(residuo, bins=30, color='#4E148C', edgecolor='black')
    ax_hist.set_title(f"Histograma de los Residuos - ventana {k}")
    ax_hist.set_xlabel("Valor del Residuo")
    ax_hist.set_ylabel("Frecuencia")

    sm.graphics.tsa.plot_acf(residuo, lags=7, ax=ax_acf, color='#4E148C', vlines_kwargs={'color': '#4E148C'}, alpha=0.05)
    ax_acf.set_title(f"Autocorrelación de los Residuos - ventana {k}")

plt.tight_layout()
plt.show()

En las gráficas anteriores podemos ver los histogramas y autocorrelaciones para los residuos del modelo, aunque parecen seguir una distribución normal, sabemos que la prueba analítica arrojó que no se cumple el supuesto de normalidad y en las autocorrelaciones se puede ver un ligero patrón lo que evidencia la correlacion que existe entre ellos.

Predicción en el conjunto de test

sns.set(style="whitegrid")

fig, axes = plt.subplots(4, 1, figsize=(12,8))

for i, (k, datos) in enumerate(Test_xgboost.items()):
    obs = datos["real"]
    predic = datos["predic"]
    yTrain = datos["trainY"]
    scsoreTest=datos['scoretest']
    ytrain=yTrain[-400:]
    ax = axes[i]    
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(obs)), y=obs, label=f"Observado", color='#8e94f2', ax=ax)
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(predic)), y=predic, label=f"Predicho", color="#f77f00", linestyle="--", ax=ax)
    sns.lineplot(x=range(len(ytrain)), y=ytrain, label="Entrenamiento", color="#4E148C", ax=ax)
    
    ax.set_xlabel("")
    ax.set_ylabel("")
    ax.set_title(f"Ventana {k}, score en el test= {round(scsoreTest,2)}")
    ax.legend()

plt.tight_layout()
plt.show()

En las gráficas anteriores se aprecian los ajustes del modelo en los diferentes conjunto de test, vemos que las predicciones van siguiendo relativamente bien pero este modelo sigue teniendo dificultades al momento de predecir los valores extremos.

Prediccion conjunto de entrenamiento

# Crear la figura
fig = px.line(data_pm10, x='fecha', y=[y, predicc_xgb], 
              labels={'value': 'PM10', 'variable': 'Leyenda'})

# Asignar nombres a las trazas y actualizar el color y el estilo
fig.data[0].name = 'pm10_observado'  # Serie observada
fig.data[1].name = 'pm10_predicho'  # Serie predicha

# Cambiar el estilo de las líneas
fig.update_traces(line=dict(width=2))  
fig.update_traces(line_color='#9d4edd', selector=dict(name='pm10_observado'))  
fig.update_traces(line_color='#3c096c', selector=dict(name='pm10_predicho')) 

# Configurar el diseño
fig.update_layout(
    title='Predicción vs Observación de PM10: XGBoost',
    xaxis_title='Fecha',
    yaxis_title='Concentración de PM10',
    legend_title='Leyenda',
    template='plotly_white',
    hovermode='x unified',
    width=1000,
    height=500,
    legend=dict(x=1, y=1, xanchor='center', yanchor='bottom')
)

# Mostrar la gráfica
fig.show()

SVR

Veamos el rendimiento del modelo SVR en el conjunto de entrenamiento y los diferentes conjuntos de test.

Hacemos la separación de los conjuntos respectivos de test.

Test_svr={}
ventanas=[7,14,21,28]
for i in ventanas:
    Test_svr[str(i)]= {'trainX': X[:-i*24], 'trainY':y[:-i*24],'xtest':X[-i*24:],'real': y[-i*24:]}

Veamos los mejores parametros que entregan el mejor score en e conjunto de entrenamiento.

pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('svr', SVR())  
])

tscv = TimeSeriesSplit(n_splits=5)

param_grid = {
    'svr__C': [0.1, 1, 10],          
    'svr__gamma': [ 0.01, 0.1, 1],  
    'svr__kernel': ['linear', 'rbf'] 
}

# Realizar GridSearchCV
grid_search = GridSearchCV(pipeline, param_grid, cv=tscv, scoring='r2')
grid_search.fit(X, y)

# Imprimir los mejores parámetros y el score de entrenamiento
print("Mejores parámetros:", grid_search.best_params_)
print("Score de entrenamiento:", grid_search.best_score_)

# Realizar predicciones con el mejor modelo
predicc_svr = grid_search.predict(X)

# Calcular los residuos
residuos_svr = y - predicc_svr

# Guardar el modelo, las predicciones y los residuos en un archivo pickle
with open('modelo_svr.pkl', 'wb') as f:
    pickle.dump({
        'modelo': grid_search.best_estimator_,
        'predicciones': predicc_svr,
        'residuos': residuos_svr
    }, f)

print("Modelo, predicciones y residuos guardados en 'modelo_svr.pkl'")

Mejores parámetros: {'svr__C': 10, 'svr__gamma': 0.1, 'svr__kernel': 'rbf'}
Score de entrenamiento: 0.5731122924905874
Modelo, predicciones y residuos guardados en 'modelo_svr.pkl'

Veamos que la comnbinación de parametros que entrega el mejor score en el entrenamiento, siendo este un \(R^2= 57%\), la combinacion de un C=10, un kernel radial y \(\gamma =0.1\). Se debe mencionar que este score es mucho mas bajo de lo que se esperaba para este modelo, se debe revisar como mejorarlo. Veamos el rendimiento de este modelo en los distintos conjuntos de test.

with open('modelo_svr.pkl', 'rb') as f:
    model_svr = pickle.load(f)

for clave, valor in Test_svr.items():

    xtrain = valor['trainX']
    ytrain = valor['trainY']
    xtest = valor['xtest']
    real = valor['real']
    
    pipeline = Pipeline([
        ('scaler', StandardScaler()),
        ('svr', SVR())  
    ])

    tscv = TimeSeriesSplit(n_splits=5)

    param_grid = {
        'svr__C': [0.1, 1, 10],          
        'svr__gamma': [ 0.01, 0.1, 1],  
        'svr__kernel': ['rbf'] 
    }
    grid_search = GridSearchCV(pipeline, param_grid, cv=tscv, scoring='r2')
    grid_search.fit(xtrain, ytrain)

    print(f'Test {clave}')
    print("Mejores parámetros:", grid_search.best_params_)
    print("Score de entrenamiento:", grid_search.best_score_)

    predicc_svr = grid_search.predict(xtest)
    
    Test_svr[clave]['predic'] = predicc_svr
    Test_svr[clave]['residuos'] = real - predicc_svr
    Test_svr[clave]['parametros']=grid_search.best_params_

# Guardar el diccionario Test_svm en un archivo pickle
with open('Test_svm1.pkl', 'wb') as f:
    pickle.dump(Test_svr, f)

print("Diccionario Test_svm guardado exitosamente en Test_svm.pkl")

Test 7
Mejores parámetros: {'svr__C': 10, 'svr__gamma': 0.1, 'svr__kernel': 'rbf'}
Score de entrenamiento: 0.572672728981732
Test 14
Mejores parámetros: {'svr__C': 10, 'svr__gamma': 0.1, 'svr__kernel': 'rbf'}
Score de entrenamiento: 0.5731795930082803
Test 21
Mejores parámetros: {'svr__C': 10, 'svr__gamma': 0.1, 'svr__kernel': 'rbf'}
Score de entrenamiento: 0.5736043579260931
Test 28
Mejores parámetros: {'svr__C': 10, 'svr__gamma': 0.1, 'svr__kernel': 'rbf'}
Score de entrenamiento: 0.5745073227399055
Diccionario Test_svm guardado exitosamente en Test_svm.pkl

Se entrenó un modelo para cada uno de los conjuntos, con el fin de tener una mejor predicción en cada conjuntos.Se puede apreciar que los scores para cada entrenamiento no dista mucho del score obtenido de manera general, puesto que fueron escogidos los mismos parametrps. Veamos ahora si, como es el score en el conjunto de test.

with open('Test_svm1.pkl', 'rb') as f:
    Test_svr_pkl = pickle.load(f)

tabla5=pd.DataFrame()
for clave, valor in Test_svr_pkl.items():
    real = valor['real']
    predic = valor['predic']
    residuos=valor['residuos']
    
    # Calcular las métricas
    mae = mean_absolute_error(real, predic)
    mse = mean_squared_error(real, predic)
    rmse = np.sqrt(mse)
    r2 = r2_score(real, predic)
    Test_svr_pkl[clave]['scoretest']=r2

    ljung_box_results = acorr_ljungbox(residuos, lags=[10], return_df=True)
    
    jb_stat, jb_pvalue = jarque_bera(residuos)
    
    tabla_temp = pd.DataFrame({
        'Ventana': [clave],
        'MAE': [mae],
        'MSE': [mse],
        'RMSE': [rmse],
        'R2': [r2],
        'Ljung-Box p-value': [ljung_box_results['lb_pvalue'].values[0]],
        'Jarque-Bera p-value': [jb_pvalue]
    })
    
    tabla5 = pd.concat([tabla5, tabla_temp], ignore_index=True)

tabla5

	Ventana	MAE	MSE	RMSE	R2	Ljung-Box p-value	Jarque-Bera p-value
0	7	6.443400	94.956103	9.744542	0.627489	0.010996	6.315256e-175
1	14	7.014769	105.840744	10.287893	0.640914	0.020017	1.720035e-114
2	21	7.131726	103.493423	10.173172	0.597030	0.001179	9.446405e-95
3	28	7.311465	112.142289	10.589726	0.550346	0.005492	2.966446e-145

El mejor score obtenido de todas las ventanas fue nuevamente la ventana de 14 días para el test, con el mejor score que hemos visto que es de 65%, notamos que los errores son bastante bajos, pero aun así no se cumplen los supuestos de independencia y normalidad en ninguno de los conjuntos.

sns.set(style="whitegrid")
fig, axes = plt.subplots(4, 2, figsize=(8, 8))  

for i, (k, datos) in enumerate(Test_svr_pkl.items()):
    residuo = datos["residuos"]
    ax_hist = axes[i, 0]  
    ax_acf = axes[i, 1]   

    ax_hist.hist(residuo, bins=30, color='#4E148C', edgecolor='black')
    ax_hist.set_title(f"Histograma de los Residuos - ventana {k}")
    ax_hist.set_xlabel("Valor del Residuo")
    ax_hist.set_ylabel("Frecuencia")

    sm.graphics.tsa.plot_acf(residuo, lags=7, ax=ax_acf, color='#4E148C', vlines_kwargs={'color': '#4E148C'}, alpha=0.05)
    ax_acf.set_title(f"Autocorrelación de los Residuos - ventana {k}")

plt.tight_layout()
plt.show()

En las gráficas anteriores podemos ver los histogramas y autocorrelaciones para los residuos del modelo, aunque parecen seguir una distribución normal, sabemos que la prueba analítica arrojó que no se cumple el supuesto de normalidad y en las autocorrelaciones se puede ver un ligero patrón lo que evidencia la correlacion que existe entre ellos.

sns.set(style="whitegrid")

fig, axes = plt.subplots(4, 1, figsize=(12,8))

for i, (k, datos) in enumerate(Test_svr_pkl.items()):
    obs = datos["real"]
    predic = datos["predic"]
    yTrain = datos["trainY"]
    scsoreTest=datos['scoretest']
    ytrain=yTrain[-400:]
    ax = axes[i]
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(obs)), y=obs, label=f"Observado", color='#8e94f2', ax=ax)
    sns.lineplot(x=range(len(ytrain), len(ytrain) + len(predic)), y=predic, label=f"Predicho", color="#f77f00", linestyle="--", ax=ax)
    sns.lineplot(x=range(len(ytrain)), y=ytrain, label="Entrenamiento", color="#4E148C", ax=ax)
    
    ax.set_xlabel("")
    ax.set_ylabel("")
    ax.set_title(f"Ventana {k}, score en el test= {round(scsoreTest,2)}")
    ax.legend()

plt.tight_layout()
plt.show()

Lo mas importante para resaltar de esta gráfica de las predicciones hechas en cada conjunto de test realizadas por el modelo de SVR es que este modelo lográ predecir mucho mejor los extremos de la variable PM10, lo que hace que este modelo sea un poco mejor a los demas modelos, aunque los scores no sean los mejors.

# Crear la figura
fig = px.line(data_pm10, x='fecha', y=[y, model_svr['predicciones']], 
              labels={'value': 'PM10', 'variable': 'Leyenda'})

# Asignar nombres a las trazas y actualizar el color y el estilo
fig.data[0].name = 'pm10_observado'  # Serie observada
fig.data[1].name = 'pm10_predicho'  # Serie predicha

# Cambiar el estilo de las líneas
fig.update_traces(line=dict(width=2))  
fig.update_traces(line_color='#9d4edd', selector=dict(name='pm10_observado'))  
fig.update_traces(line_color='#3c096c', selector=dict(name='pm10_predicho')) 

# Configurar el diseño
fig.update_layout(
    title='Predicción vs Observación de PM10: SVR',
    xaxis_title='Fecha',
    yaxis_title='Concentración de PM10',
    legend_title='Leyenda',
    template='plotly_white',
    hovermode='x unified',
    width=1000,
    height=500,
    legend=dict(x=1, y=1, xanchor='center', yanchor='bottom')
)

# Mostrar la gráfica
fig.show()